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Abstract 



D ■ This paper presents a fortran program to solve diverse few-body problems 

" . with the stochastic variational method. Depending on the available compu- 

tational resources the program is applicable for A^ = 2 — 3 — 4 — 5 — 6 — ...- 
body systems with L = total orbital momentum. The solution with the 
stochastic variational method is "automatic" and universal. One defines the 
^ ■ system (number of particles, masses, symmetry, interaction, etc. ) and the 

cn ■ program finds the ground state energy and wave function. The examples in- 

iH , elude nuclear (alpha particle: four-body, ^He: six-body), atomic (td^^ and 

f^ I e'^e~e~^e~) and subnuclear (the nucleon and the delta in a nonrelativistic 

^ I quark model) systems. The solutions are accurate for excited states as well, 

and even the Efimov-states can be studied. The program is available from 
the author (e-mail:varga@rikaxp. riken.go.jp). 

T3 ; PACS number(s): 27.20.-Fn, 23.20.-g, 23.40.-s, 21.60.Gx 

^ ' 



Typeset using REVT^ 



I. INTRODUCTION 

The solution of few-body problems is an important basic problem of physics. One en- 
counters few-body problems from celestial to quark level in atomic, chemical, nuclear or 
subnuclear physics. In recent years, due to the intense experimental, theoretical and tech- 
nological interest in mesoscopic scale systems in solid state physics (few ions in a trap, few 
electrons in a quantum dot, etc.) the traditional domain of application is greatly enlarged. 

In the last few years we have elaborated a powerful technique, the stochastic varia- 
tional method (SVM) [0,|],|ll, which is proved to be especially suitable for solution of diverse 
few-body problems. The stochastic variational method optimizes the variational basis in 
a random trial and error procedure. The basis selection is free from any bias, keeps the 
dimension of the basis low, and lastly but most importantly, provides a very accurate solu- 
tion. The method can be used with different type of bases. The correlated Gaussian basis 
^ seems to be particularly suitable for description of A^ = 2 — 8 particle systems, while for 
larger number of particles displaced Gaussian or harmonic oscillator bases can be applied. 

The method is a natural extension of the rigorous few-body (A^ = 3, 4) techniques [5-11] 
to larger systems of strongly correlating particles and offers a wide range of applications. 

The aim of this paper is to present a computer code for solution of few-body problems 
with the stochastic variational method on correlated Gaussian bases. The particles can 
interact via different central (Coulomb, Yukawa, Gauss or power law or other numerically 
given) potentials. The interaction may contain spin-isospin dependent operators. The pro- 
gram is general: the number of particles is in principle arbitrary, in practice it is limited 
by the speed and the memory of the available computer. The method can be used to get 
very accurate solution for smaller systems or to find an approximate upper bound for larger 
systems. One can treat bosons and fermions. 

The applicability of the program is shown with various examples, including nuclear, 
Coulombic and quark systems. The accuracy of the solutions are tested by comparing the 
results to those of the literature. An example of the Efimov-states |T^ shows that the 
method gives precise energies for ground and excited states. 

The plan of the paper is as follows. In section 2 we outline the method. In section 3 
we show the calculation of the matrix elements. The fortran code is described in section 4. 
Examples are presented in section 5. The papers ends with a brief summary. 

II. THE STOCHASTIC VARIATIONAL METHOD 

Let's consider an A^-particle system, where the ith particle with mass rrii, spin Sj, isospin 
ti and charge Zt is placed at the position Tj. The positions of the particles can be more 
conveniently described by introducing a set of relative (Jacobi) coordinates x = (xi, ..., x^r.i) 
and the coordinate of the centre of mass xat. The object of this paper is to solve the many- 
body Schrodinger-equation 

iV 2 N 

for central two-body interactions Vij. 



In our variational approach the basis functions are assumed to have the form 

IpSMsTMri^, ^) = '^{G'a(x)X5Ms?7tMt}, G'a(x) = e^^''^'', (2) 

where the operator A is an antisymmetrizer, Xsms is the spin function, and tjtmt is the 
isospin function of the system (for Coulombic system this latter can be suppressed). The 
diagonal elements of the {N — 1) x {N — 1) dimensional symmetric, positive definite matrix 
A corresponds to the nonlinear parameters of an Gaussian expansion, and the off diagonal 
elements connect the different relative coordinates representing the correlations between 
the particles. This trial function, the correlated Gaussian basis, is widely used in physics 
H JTSlJT^ , although the applications are mostly restricted to three and four particle systems. 
The above form assumes orbital angular momentum L = 0. 

In the variational method the wave function of the system is expanded as 

K 

^ = ^ Cil/jsMsTMr (x, Ai), (3) 

i=l 

and an upper bound for the ground state energy of the system is given by the lowest eigen- 
value of the generalized eigenvalue problem 

HC = EkBC, (4) 

where 

Hij = ii'SMsTMr (x, Ai) \H\iJsMsTMT (x, Aj)) , Bij = {ipSMsTMr (x, Ai) HjsMsTMt (x, Aj)) 

(5) 

The adequate choice of the nonlinear parameters (the elements of the {N — 1) x {N — 1) 
matrix A) is crucial. If these parameters are properly optimized, the variational solution 
gives very precise energies. Although several ways for the choice of elements of A have been 
suggested, there is no safe recipe available. While the numerical optimization would be in 
principle the method of choice, in practice there are several difficulties to face. The most 
serious ones amongst these are probably the large number of parameters to be optimized and 
the nonorthogonality of the basis states. Due to the nonorthogonality, none of the parameter 
sets is indispensable, and several different choices can represent the wave functions equally 
well. This property makes the optimization tedious but offers the possibility of random 
selection of the nonlinear parameters. 

In the stochastic variational method the most suitable parameters are selected as follows: 

(1) several sets of (A", ....,A^,n = 1, ...,M) are generated randomly, 

(2) by solving the eigenvalue problem the corresponding energies (-E^, ■■■E^) are determined, 

(3) the parameter set (A", ...., A^) belonging to the lowest energy E]^ are chosen as basis 
parameters. 

The bottleneck of the above procedure is that it assumes a full diagonalization to solve the 
eigenvalue problem and beyond a certain dimension it becomes too computer time consum- 
ing. In practical implementation it is more advantageous to fix (Ai, ....,Ak-i) and change 
only A^. In this case only one row (column) of the matrices H and B has been changed, 



and no diagonalization is needed but one can utilize the formulae of ref. [3]. Moreover, as 

by increasing the number of basis states the energy decreases (-Ex+i < E^), we found it 

advantageous to increase the dimension of the basis after selecting the most suitable A]^ 

amongst the random candidates. This step ensures that the energy E^, belonging to a new, 

randomly selected (A^) basis state, is lower than the energy Ek~i on the previously found 

basis, and provides a convenient selection criteria. The steps of the most economical way is 

therefore: 

(i) several sets of (A^,n = 1, ...,A/') are generated randomly, 

(ii) by solving the eigenvalue problem the corresponding energies {E}^, ■■■E^) are determined, 

(iii) the parameter {A^) belonging to the lowest energy E^ chosen as a basis parameter 

and by adding it to the previous basis states the parameters of the new basis become 

{Ai,....,Ak-uAk = A'I^), 

(iv) the basis dimension is increased to -ft' + 1. 



III. MATRIX ELEMENTS 

One of the main advantages of the correlated gaussian basis is that its matrix elements 
can be easily calculated analytically. We list these matrix elements in the following. 
The overlap of the correlated Gaussians takes the simple form: 



{Ga\Gai) 



(27r 
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^det{A + A)) 
The matrix element of the kinetic energy between the correlated Gaussians reads as: 



(6) 
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Tcm\GA') = {Ga\Ga') 3Tr(AA) - 3Tr(A + A)- {A'AA 



(7) 



where A is an (N — 1) x (A^ — 1) diagonal matrix 
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and the reduced masses are given by 
mi+imi2-i 



/ij 



t = l A^-1 



mu-i+i 



and /ijv = 1TI12...N, 



(9) 



with mi2...i = mi+m2+ ■ ■ ■ +mi. 

To avoid the dependence of the matrix elements of the two-body interaction on the 
specific form of the potential, it is advantageous to express the potential in the form 



V{ri - Tj) = / drV{r)6{ri - r^ 



r . 



(10) 



To calculate the matrix elements of the potential we first express the relative distance vector 
Fj — Tj by the Jacobi coordinates. The Jacobi coordinates are defined as 
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The relative distance vector between two particles then can be written as 



N-l 



^i -^j= J2 Bijk^k, Bijk = Uij} - Uj)}. 



i=l 



By using this expression the matrix element of the potential is given by: 



{GA\S{r,-r,-r)\GA') = {GA\GA') 



2npi 



g 2pjj 



where 



k=l 1=1 



(11) 



(12) 



(13) 



(14) 



(15) 



By integrating eq. (14) over r one can recover the norm of the wave function. 
To calculate the matrix element of the potential one has to multiply eq. (14) by the radial 
form V"(r) and integrate over r: 



{Ga\V^,\Ga') = J drV{r){GA\S{r,-r^ -r)\GA') = {G a\G A')v{P^J) , 



(16) 



where 



v{p 



ijj 
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drV{r)e 



2pij 



(17) 



In the applications the two-body interaction is expressed as 



No / Nt \ 

V^J = J:[J:Vki\r^-r,\)]0, 
j=i \fc=i / 



(18) 



where the operators {Oi, . . . , O4) are the Wigner, Majorana, Bartlett and Heisenberg oper- 
ators. 



If the radial part can be given in the form 

V^{r) = r"e""'''+^"(n > -2), (19) 

then the integration can be calculated analytically by using the formula 

- ,„..,-.... ^ 1 ,_1). g ^-^fimn - k), (20) 

/ 1 \ *^ [fc/2] 7,1 

g{0) = eTfc{y), g{k) = {-1)'^ (-^] H,^M, {k>l), y = ^. (22) 

In other cases one has to rely on a numerical integration. To calculate the matrix elements, 
the function v{p) has to be evaluated many times for many different values of p. Both the 
analytical and the numerical evaluation takes some time on the computer. This part of 
the computation can be made faster, however, noticing that v{p) is a rather simple smooth 
function of p and it can be easily interpolated. To this end, for a given potential we tabulate 
v{p) at certain representative values of p and during the computation of the matrix elements 
v{p) is interpolated in the necessary points. The precision of the interpolation can be easily 
controlled. 



IV. SYMMETRIZATION 

The antisymmetrizer A is defined as 

N\ 

A = J2p,P^, (23) 



i=l 



where the operator Vi changes the particle indices according to the permutation {p\, ...p]v) 
of the numbers (1, 2, ..., A^), and pi is the parity of that permutation. The effect of this 
operator on the set of position vectors (ri, ..., r^) is 

Pi(ri,...,r^) = (rp.,...,rp>^) (24) 

By representing the permutations by the matrix 

(Ci)^^. = 1 if j = pi and (d),^. = otherwise, (25) 

(for example, if iV = 3 the permutation (3 1 2) is represented by 



C = 1 , (26) 




while for (1 2 3) C is a unit matrix), the effect of the permutation operator on the single 
particle coordinates 

P,(ri,...,r^)=a(ri,...,r^). (27) 

By using eqs. (11) and (27) the permutation of the relative coordinates is expressible as 

P.x = P,x, (28) 

where Pi is an {N — 1) x [N — 1) matrix obtained by omitting the last row and column 
(corresponding to the permutation invariant center-of-mass coordinate) of the N x N matrix 
U^QU. 

The correlated Gaussian function, after permutation takes the form: 

V^Ga{^)=Gpt^pX^). (29) 

This simple transformation property under permutation is particularly useful in calculating 
the matrix elements of the antisymmetized basis functions. 

In the spin-isospin space the permutation operator interchanges the indeces the single 
particle spin-isospin functions and can be easily evaluated. 

V. PROGRAM DESCRIPTION 

A. Fortran code 

To prepare the data for the calculation, the main program calls the following subroutines: 
input_data: This subroutine reads the input files described in the previous subsection. 
cm_rel_tr: This subroutine is responsible for the separation of the relative and centre-of- 
mass coordinates. 

me_spiso: The overlap of the spin-isospin part, and the Pi matrices are calculated in this 
segment. 

The matrix elements of the kinetic energy and potential energy are calculated in the 
subroutines vkin_ene and vpot_ene by using the formulae (7) and (16). The overlap of 
the basis states, given by eq. (6), is computed in vove_mat. The subroutine inat_elem 
prepares the matrices H and B of eq. (4) . 

The parameters of the basis states can be selected in three ways. The first (si) is a direct 
calculation on a given basis, the second (s2) and third (s3) are two versions of the stochastic 
selection of the basis states. 

(si) One can use predefined basis. In this case, the parameters of the basis has to be written 
in the file "fbs.res". The first line of this file defines the number of basis states K. In 
the next lines the parameters of the basis are listed. Each line contains the corresponding 
nonlinear parameters {Ak)ij {k = 1, ...,K). 

(s2) One can select the parameters randomly and set up the basis step by step increasing the 
dimension, as described by (i)-(ii)-(iii)-(iv). This is a rather automatic procedure, the only 



thing one has to define before starting is the interval from which the uniformly distributed 
random numbers are generated. 

(s3) One can select the basis states through the steps (i)-(ii)-(iii), on a fixed basis dimension 
K. Before starting this procedure, one has to give initial values of the nonlinear parameters. 
Then a basis state, say the last one (the i^th) is subjected to the procedure (i)-(ii)-(iii). If a 
better new parameter set is found, then the Kih state is substituted by the new parameters. 
This procedure is repeated several times cyclically for the other elements of the basis as 
well. We refer to this step as refinement because in this way we can further improve a 
basis. Particularly, the bases, selected in the (si) or the (s2) manner can be refined with 
this method. 

The subroutines preset, svml, svm2 corresponds to the three possibilities above. To 
select amongst the three possibilities the value of ico should be set to 1,2 or 3. 

One does not need to solve the whole eigenvalue problem in each time when a basis 
parameter is changed. The program utilizes the simplifications described in ref. [3]. 

The function "f (p)" in eq. (17) can be calculated in three ways. 

(pi) One can numerically integrate over the radial coordinate in eq. (17) for given values of 

p and tabulate the function v{p). This tabulation can be done before the real calculation is 

started. During the basis selection, one needs the value of v{p) at arbitrary values of p and 

these values will be approximated by interpolation. 

(p2) If the potential is linear combination of terms like in eq. (18) one can use the analytical 

formula (19). 

(p3) In this case we numerically integrate over the radial variable in eq. (17) for each value 

of p that appear during the basis selection. The difference between this case and the first 

one is that as the actual value of p is available only during the basis selection process, this 

numerical integration can only be done then. 

The order of the above three possibilities refiects their speeds. The first way is the fastest 
by far, but it accumulates the inaccuracies of the numerical integration and interpolation. 
In most of the practical cases these inaccuracies can be kept under control. The second 
case is slower but it is exact for Coulomb, linear, harmonic oscillator, Gaussian and Yukawa 
potentials. The third case is quite slow. The recommended way is to use the first possibility 
for the basis selection. Once the selection is finished, the chosen basis can be considered as 
a predefined basis and one can rerun the calculation on that basis recalculating the matrix 
elements by the (pi) or (p2) method. The accuracy of the numerical integration is controlled 
by the variable eps. To check the accuracy of the numerical integration one may rerun the 
calculation for different values of eps. 

To select amongst (pl)-(p2)-(p3) one has to create a file "pot.inp" and write a number 
"ipcon" (ipcon=l,2 or 3) into its first line. In the case (p2) one should write the parameters 
of the potentials into the file "pot.inp" as well. The number of the linear combinations 
of terms like eq. (17) should be written into the second line. It is to be followed by the 
coefficient of the linear combination, and the parameters a, 6, n of the terms in the successive 
lines, in turn. In the cases (pi) and (p3), the radial form of the potential should be defined 
in the subroutine "pot" by the user. The value of "pot" should be equal to the value of the 
two-body potential at the radial distance "x" . 



B. Selection of nonlinear parameters 

The Jacobi coordinates are convenient, because the kinetic energy can be expressed 
simply (eq. (1)). The choice of the nonhnear parameters, however, is simpler in a coordinate 
system, where the interparticle distances Tj — r^) are used. The basis function in that system 
takes the form: 

exp{-5Iaii(ri-rj)^}. (30) 

The elements of A in G'^(x) can be expressed by aij using an appropriate linear trans- 
formation. The advantage of a-ij is that it is more directly connected to the interparticle 
distances and it can be more uniformly used. In the practical applications Uij are generated 
as random numbers from the \bmini bmax] interval. The elements of A are then expressed by 
using the values of aij. 

Although one can choose the elements of A (through aij) independently of each other 
randomly, we have found that it is slightly more advantageous to follow the following way. 
At first, one generates random values for each element. Then the first element of the matrix 
changed randomly Kq times. The best (giving the lowest ground state energy) parameter 
is selected and fixed. The second element is fixed in the same way and this process is 
continued till the last element. This procedure is then repeated Mq times. This selection 
requires A'"o x Mq evaluations. Only one element of the matrix of the nonlinear parameters 
is changed in each step and that leads to the possibility of fast evaluation. 



C. Input data 

The input data specifies the system and the interaction. One has to define the number of 
particles, the spins, isospins, charges and masses. The spin and isospin states are represented 
by integers, the 'down' state is coded as 1 and the 'up' state is coded as 2. The input data 
is read from the input file "fbs.inp". The first line defines the number of particles A^, 
the second line lists the masses of particles {mi, . . . ,m]\f), the third line gives the charges 
(zi, . . . , zjy), the fourth the isospins (ti, . . . , tjy). The spin part of the system can be given as 
linear combination of spin states. The next line contains the number of linear combinations, 
and the succeeding lines give the linear combination coefficients and the spin states. The 
following line should contain the value of h /m, an initial (negative integer) number for the 
random number generator and the control parameters ico and ibf. The first parameter selects 
the method of basis optimization (sl),(s2) or (s3), while the second defines the symmetry 
of identical particles (ibf=l for fermions and ibf=2 for bosons). The next line gives Mq, Kq 
and K. The program generates random numbers for a given nonlinear parameter, say An 
Kq times, and once the best parameters for all Aij has been found this random selection is 
repeated Kq times. 

In the fortran program, the parameter "mnbas" and "mnpar" defines the maximum 
number of basis states and the maximum number of particles, respectively. 



D. Output data 

The main results of the calculation are the energies of the ground and excited states and 
the wave function. The energies are written in the file "ener.dat". The basis dimension, the 
coefficients of the wave function of the ground state and the nonlinear parameters of the 
basis are in the file "fbs.res" after the calculation is finished. 



VI. EXAMPLES 

In this part we show a few examples for the application of the program. The compu- 
tational time varies for different cases and strongly depends on the basis size and on the 
number of particles. 

A. The fd/U and the positronium molecule 

This example includes a three- and a four-body system with Coulomb interaction. The 
tdfi" molecule is a bound system of two positively and a negatively charged particles with 
unequal masses. This system attracted considerable attention in connection with the muon 



catalysed fusion ||T5[. The input parameters (the input files) are shown in Table I. Atomic 
units are used. The isospin quantum number is not necessary in these systems so it just 
stands for distinguishing of the particles with different charges. The matrix elements of the 
Coulomb-potential can be analytically calculated and thus the second way (p2) is used. 

One can start the calculation by a step-by-step random selection of the basis states (s2- 
type, ico=2). The first four digits of the ground state energy can be reproduced on a basis 
size of i^ = 50 (see Table II.). To improve the energy on this basis size further, one is to 
restart the calculation with a refinement circle (s3-type, ico=3). By repeating the refinement 
three times, the fifth digits of accuracy can be reached. By increasing the basis size more 
accurate result can be found in need. The energies in the tables are written into the output 
file "ener.dat". 

One can solve the two-body problem as well. The above example can be recalculated 
for two-body case (t/i~), by changing the number of particles to 2 in the input file. One 
immediately arrives at the energy of the (t/i~) threshold (-99.64 a.u.). 

The next example is the positronium molecule, the e~^e~e~^e~ Coulomb four-body system. 
This system has been subject of quite a few study, but accurate calculations have become 
available only recently PJT^JT^ . In this case we have four particles with equal masses and a 



nonadiabatic treatment is necessary. The spins of the particles are coupled to 5* = 0. This 
spin coupling leads to a more symmetric system and the convergency is much faster. Due 
to the mass difference between the particles in this and in the previous case, the interval of 
the random numbers are to be chosen differently. The input files are in Table III. The basis 
selection follows the same way as in the previous example and the results are compared in 
Table IV. to those of other calculations. To increase the accuracy further one has to use 
a larger basis. It is remarkable that this calculation, after refining cycles, gives very close 
results on a basis dimension of if = 50 to that of refs. [1,13,14] with K = 300 Correlated 
Gaussians. 

10 



To see the tremendous effect of the spin couphng on the convergency one may try to 
rerun the same calculation with using only one of the spin configurations. 

B. The alpha particle and the ^He 

In this subsection we describe the application of the program for nuclear systems with a 
simple central interaction. A four- and a six-body system (the alpha particle and ^He) are 
chosen as example. 

In the first example of this subsection the ground state energy of the alpha particle is 
calculated with the central, spin-isospin independent Malfieit-Tjon V |T^ (MTV) interac- 
tion. The input is given in Table V. The spin coupling is the same as in the case of the 
positronium molecule. Although the potential matrix elements are analytical as before, the 
first representation (pi) is suggested as an illustrative example. To this end one has to write 
a function (in fortran) "pot(r)" which defines the MTV potential. The program runs faster 
with this choice and the accuracy is still appropriate for nuclear systems. The results are 
compared to other calculations in Table VI. 

The next example is a six-body system. Usually only the Quantum Monte Carlo [5,10] 
methods are capable of going beyond four particles. In this example we try to solve the ^He 
with the central, spin-isospin dependent Minnesota potential |T^. The potential, though 
again analytical, is used in numerical form as a fortran function. To illustrate another option, 
this example uses a predefined basis (ico=l, sl-type). The predefined basis is to be included 
in the file "fbs.res". (The program package contains this file under the name "fbsl.res".) 

C. The nucleon and the delta in a nonrelativistic quark model 

In this example the SVM is used to solve the nonrelativistic three-quark problem. The 
nonrelativistic three-quark model of baryons has a long history (see e.g. refs. [18-21]). Many 
of the properties of the baryons are quite successfully explained in this framework ||T8| , p!9| , ^ . 
In this section we solve the three-body Schrodinger with a one-gluon exchange potential ||2D 



for the nucleon and the delta. The parameters of the potential is taken from ref . [20] . The 
three-quark problem with this potential has been studied in the framework of the Faddeev- 
equations in ref. [21]. We use the same parameters and compare our results to theirs. 

The potential has no dependence on color degrees of freedom therefore the antisymmetric 
color part of the wave function can be factorized and the spin-flavor-space part of the wave 
function must be symmetric (the control parameter "ibf" should be set to 2). The input 
files for the nucleon and the delta differ in the spin- and isospin- (flavor) part (see Table ). 
The potential is given in a numerical form (ipcon=l). The quark masses is taken to be 337 
MeV like in ref. [21]. 

Due to the confining potential the basis size required for a reasonably accurate solution 
is very small. After suitable modifications of the symmetrization of the wave function the 
program can be used to calculate other baryons as well [^ . 

11 



D. The Efimov states 



The last example shows the application for Efimov states |T^. A three-body system 
with short range forces may have several bound states. If the two-body subsystem have 
just one bound state whose energy is (close to) zero, the three-particle can interact at long 
distances and an infinite number of bound state may appear in the three-body system. These 
"Efimov states" are extremely interesting from both experimental and theoretical points of 
view because of their distinct properties. These states are very loosely bound and their 
wave functions extend far beyond those of normal states. By increasing the strength of the 
two-body interactions these states disappear. 

A three-boson system is considered, for simplicity. In this case we have three spinless 
identical particles with a symmetrized spatial wave function. The potential between the 
particles is taken as a the Posch- Teller interaction because this potential is analytically 
solvable for the two-body case, and therefore the two-body binding energy can be easily set 
to zero. A spatially very extended basis is needed to represent the weakly bound states. 
In this case we use a predefined basis (si). The predefined basis is created by using only 
the diagonal elements of the matrix A. These diagonal elements are taken as geometric 
progression: 

i^k)ii = -, ^ny-; {Ak)22 = - ^zT)-^ (i = l,...,n, j = l,...,n, k = {n-l)i + j). 



[aoqo ) loogo 



(32) 



This construction defines an n x n dimensional basis. The parameters are chosen as n = 
20, oq = 0.1, Xq = 2.4. (This predefined basis can be found in the file "fbs2.res".) To try this 
example one has to use the potential function "pot" with the Posch- Teller potential ||23|| : 



_ 625.972 1251.943 

^^^ ~ sinh(1.586 r)^ ~ cosh(1.586 r)2' ^ '' 

The result of this calculation is given in Table XL One can find the first three Efimov states 
with this basis size. The ratios of the energies of the bound states follow the Ei+i/Ei = 
exp{— 27r} rule [|12|. By increasing the basis size one can reveal more bound excited states. 



Note that the value of OoXq roughly corresponds to the spatial extension of the basis. The 
present basis in this example goes out up to oq Xq = 0.1 x 2.4^ = 1674990 (fm). This 
extension is necessary to get the third bound state. As a comparison, to calculate the 
ground state energy (—4.81) it is enough to choose n = 7 and the basis covers only the 
[0, 10] (fm) interval. The root mean square radii of the ground state is about 1.5 fm, while 
that of the first excited state is about 6000 fm. These facts illustrate the tremendous spatial 
extension of the Efimov states. 



VII. SUMMARY 

A fortran program is presented which solves the few-particle Schrodinger- equation by 
using the stochastic variational method. The usefulness and applicability of the method is 

12 



illustrated on various examples. The program can be used to solve diverse few-body prob- 
lems. The stochastic variational method selects the most important basis states and keeps 
the basis size low even for six-body problems. The solution with the SVM is "automatic" 
and universal. One defines the system (number of particles, masses, symmetry, interaction, 
etc. ) and the program finds the ground state energy and wave function. The refining steps 
[24,25] applied here increase the accuracy further. The number of particles can be easily 
increased up to a certain (computer dependent) limit. The program can be used without 
the SVM, by using a predefined basis as well. 

This work was supported by OTKA grant No. T17298 (Hungary) and by Grants-in-Aid 
for Scientific Research (No. 05243102 and No. 06640381) and for International Scientific 
Research (Joint Research) (No. 08044065) of the Ministry of Education, Science and Culture 
(Japan). K. V. gratefully acknowledges the support of the JSPS. We are thankful for the 
possibility of using the computer facilities of RIKEN. 
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TABLES 
TABLE I. Input files for tlie tdfi system. 



fbs.inp 



3 number of particles 

206.786,5496.918,3670.481 masses of particles 

—1,1,1 charges of particles 

1,2,2 isospins 

1 number of spin configurations 

1., 1,2,1 coefficient, spins 

1.0,-14491,2,1 n^/m,iran, ico, ibf 

5,100,50 Mo,Ko,K 

0.0000001,0.1 b„,in,bmax 



pot.inp 



2 ipcon 

1,1 No,Nt 

1,0.,0.,— 1 parameters of the potential 



TABLE II. Results for the tdn system. (In atomic units.) 

Energy basis selection 

-111.291 a.u. {K = 50) (s2)-type 

-111.342 a.u. {K = 50) (s2)-type followed by an (s3)-type 

—111.360 a.u. {K = 50) (s2)-type followed by 3 times (s3)-type 

-111.36444 a.u. {K = 200) ref. [25] 

-111.36451 a.u. {K = 1442) ref. [15] 
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TABLE III. Input files for the positronium molecule. 



fbs.inp 



4 number of particles 

1 . , 1 . , 1 . , 1 . masses of particles 

1,1,-1,-1 charges of particles 

1,1,2,2 isospins 

4 number of spin configurations 

1., 1,2, 1,2 coefficient, spins 

-l.,l,2,2,l 

-l.,2,l,l,2 

l.,2,l,2,l 

1.0,-14491,2 ti^/m^ran, ico, ibf 

5,25,50 Mo,Ko,K 

U.i,iO. OjYiintOuiax 



pot.inp 



2 ipcon 

1,1 No,Nt 

1,0.,0., — 1 parameters of the potential 



TABLE IV. Results for the positronium molecule. (Atomic units are used.) 

Energy basis selection 

-0.51548 a.u. {K = 50) (s2)-type 

—0.51575 a.u. (K = 50) (s2)-type followed by an (s3)-type 

—0.51586 a.u. (K = 50) (s2)-type followed by 3 times (s3)-type 

-0.51600 a.u. (K = 300) refs. [1,13,14] 
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TABLE V. Input files for the alplia particle. 



fbs.inp 

4 number of particles 

1 . , 1 . , 1 . , 1 . masses of particles 

1.1.1.1 charges of particles 

1.1.2.2 isospins 

4 number of spin configurations 

1., 1,2, 1,2 coefficient, spins 

-l.,l,2,2,l 

-l.,2,l,l,2 

l.,2,l,2,l 

41.47,-14491,2 h^/m^ran, ico, ibf 

5,25,50 Mo,Ko,K 



0.1,15. 



-'mini "max 



pot.inp 



1 ipcon 

1,1 No,Nt 

0.,0.,0.,0 parameters of the potential 



TABLE VL Results for the alpha particle. (In MeV.) 

Energy basis selection 

-31.30 {K = 100) (s2)-type 

-31.32 {K = 100) (s2)-type followed by an (s3)-type 

-31.33 {K = 100) (s2)-type followed by 3 times (s3)-type 

-31.36 (K = 150) ref. [3] 
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TABLE VII. Input files for the ^He. 



fbs.inp 



number of particles 
masses of particles 
charges of particles 
isospins 

number of spin configurations 
coefficient, spins 
fi^/m^iran, ico, ibf 
Mo,Ko,K 



6 

1,1,1,1,1,1, 

1,1,2,2,2,2 

1 

l.,l,2,l,2,l,2 

41.47,-14491,2 

5,25,50 



0.1,15. 


















"mini bmax 




pot.inp 


1 

1,1 
0.,0.,0.,' 



















ipcon 

parameters 


of the potential 








TABLE VIII. 


Results 


for the 


6He. 


(In 


MeV.) 




Energy 




















basis selection 


-30.00 
-30.07 


{K-- 


= 100) 
= 300) 
















(sl)-type 
ref. [3] 






TABLE IX. 


Input files for the nucleon and the delta particle 




fbs.inp 


(nucleon) 



















l.,l.,l. 
1,1,1 

1,2,2 

2 

l.,l,2,l 

-1., 1,1,2 

115.54,-14491,2,2 

5,25,50 



number of particles 

masses of particles 

charges of particles 

isospins 

number of spin configurations 

coefficient, spins 

h /m,iran, ico, ibf 
Mo,Ko,K 



0.1,8. 


"mini bmax 


fbs.inp (delta, spin-isospin part only 


1,1,1 


isospins 


1 


number of spin configurations 


l.,l,l,l 


coefficient, spins 





TABLE X. 


Results 


for nucleon and delta 


(in MeV) (3*337 MeV 


is added). 


Energy 












basis selection 


Nucleon 


1021 {K = 
1024 


= 10) 










(s2)-type 
ref. [21] 


delta 


1330 {K = 
1330 


= 5) 










(s2)-type 
ref. [21] 






TABLE XL 


Input files for the Efimov states. 




fbs.inp (nucleon) 



3 

l.,l.,l. 

1,1,1 

1,1,1 

1 

l.,l,l,l 

41.47,-14491,1,2 

5,25,50 



dimension 
energy (1) 
energy (2) 
energy(3) 



number of particles 

masses of particles 

charges of particles 

isospins 

number of spin configurations 

coefficient, spins 

ti /m,iran, ico, ibf 

Mo,Ko,K 



0.1,15. 


"mini Umax 


pot.inp 


1 

1,1 
0.,0.,0.,0 


ipcon 

No,Nt 

parameters of the potential 


output: ener.dat 



400 

-4.81 

-8.77 X IQ-^ 
-1.50 X 10"^ 
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